TESTING HIGHER-ORDER LAGRANGIAN 
PERTURBATION THEORY 
AGAINST NUMERICAL SIMULATIONS - 



2. HIERARCHICAL MODELS 

by 

A.L. Melott 1 , T. Buchert 2 , A.G. Weifi 2 

department of Physics and Astronomy 
University of Kansas 
Lawrence, Kansas 66045 
U. S. A. 

2 Max-Planck-Institut fiir Astrophysik 
Postfach 1523 
85740 Garching, Munich 
F. R. G. 

submitted to Astronomy and Astrophysics, Main Journal 







Testing higher-order Lagrangian perturbation theory 



Summary: We present results showing an improvement of the accuracy of perturba- 
tion theory as applied to cosmological structure formation for a useful range of scales. 
The Lagrangian theory of gravitational instability of Friedmann-Lemaitre cosmogonies 
investigated and solved up to the third order in the series of papers by Buchert (1989, 
1992, 1993), Buchert & Ehlers (1993), Buchert (1994), Ehlers & Buchert (1994), is 
compared with numerical simulations. 

In this paper we study the dynamics of hierarchical models as a second step. In the first 
step (Buchert, Melott and Weifi 1994) we analyzed the performance of the Lagrangian 
schemes for pancake models, i.e., models which initially have a truncated power spec- 
trum. We here explore whether the results found for pancake models carry over to 
hierarchical models which are evolved deeply into the non-linear regime. We smooth 
the initial data by using a variety of filter types and filter scales in order to determine 
the optimal performance of the analytical models, as has been done for the "ZePdovich- 
approximation" - hereafter TZA - (as a subclass of the irrotational Lagrangian first- 
order solution) in previous work (Melott et al. 1994a). We study cross-correlation 
statistics employed in previous work for power-law spectra having indices in the range 



We find that for spectra with negative power-index the second-order scheme performs 
considerably better than TZA in terms of statistics which probe the dynamics, and 
slightly better in terms of low-order statistics like the power-spectrum. In cases with 
much small-scale power the gain from the higher-order schemes is small, but still mea- 
surable. However, in contrast to the results found for pancake models, where the higher- 
order schemes get worse than TZA at late non-linear stages and on small scales, we here 
find that the second-order model is as robust as TZA, retaining the improvement at 
later stages and on smaller scales. In view of these results we expect that the second- 
order truncated Lagrangian model is especially useful for the modelling of standard dark 
matter models such as Hot-, Cold-, and Mixed-Dark-Matter. 
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1. Introduction 



The non-linear modelling of inhomogeneities initiated by small-amplitude fluctuations at 
the time of decoupling of matter and radiation has been, for a long time, the privilege 
of numerical N-body simulations. Recently, the interest in the question has grown how 
well this modelling can be achieved by analytical or semi-numerical methods. Most of 
those efforts concentrate on Lagrangian methods introduced in a pioneering work on the 
Zel'dovich (1970, 1973) approximation (hereafter ZA). Here, we pursue a systematic ana- 
lysis of analytic solutions obtained in the framework of the Lagrangian instability picture 
investigated by one of us (Buchert 1989, 1992). Lagrangian perturbation solutions are 
available up to the third order (Buchert 1994). In this paper we restrict the presentation 
of results to the first and second-order schemes (Buchert & Ehlers 1993, Buchert 1993a). 

Parallel efforts in this direction, though less systematic, have been undertaken by 
Moutarde et al. (1991), Bouchet et al. (1992), Gramann (1993), and Lachieze-Rey 
(1993a,b). Other approximations are also built in a similar spirit, for instance the "adhe- 
sion model" (Gurbatov et al. 1989) introduces an artificial viscosity term into Zel'dovich's 
formalism, the "frozen-flow-approximation" (Matarrese et al. 1992) assumes constancy 
of the velocity field, the "frozen-potential-approximation" (Brainerd et al. 1994, Bagla 
& Padmanabhan 1994) assumes constancy of the gravitational potential. All of these ap- 
proaches are able to some extent to mimic the results of N-body simulations (for the first 
see Kofman et al. 1992 and Melott et al. 1994c, for the second, Melott et al. 1994b, for the 
third see either reference) . All of them attempt to repair a major shortcoming of Z A which 
predicts uncomfortably thick pancakes after shell-crossing: In the "adhesion model" a vis- 
cosity term mimics a 'sticking' of pancakes. In the "frozen-flow-approximation" only the 
first-generation pancakes are reached and only asymptotically. The "frozen-potential- 
approximation" (as a semi-numerical model) counteracts the thickening of pancakes by 
modelling secondary shell-crossings, a feature which is observed in numerical simulations 
as well as in the higher-order Lagrangian schemes analyzed here (compare Buchert & 
Ehlers 1993), and can therefore be attributed to part of the "true" effect of self-gravity. 
However, it is almost as time consuming as a full (Particle-Mesh) N-body simulation. 
Other work to be mentioned in this context are the articles by Matarrese et al. (1993), 
Kasai (1993), Croudace et al. (1994), Bertschinger & Jain (1994), Berschinger & Hamilton 
(1994), Kofman & Pogosyan (1994), Salopek et al. (1994), dealing with the corresponding 
equations of General Relativity, or their Newtonian limit, respectively, and Giavalisco et 
al. (1993), Munshi & Starobinsky (1994), Munshi et al. (1994). Matarrese et al. (1994) 
presented a general relativistic second-order Lagrangian solution. 
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There are two major problems with most of the tests (if there were any tests) of 
various schemes discussed here: (1) Some of them were tested by different techniques for 
accuracy, and against different simulations. This makes it impossible to compare their 
accuracy. (2) Often CDM spectra were used (e.g., Bagla & Padmanabhan 1994). Since 
the spectral index of CDM is negative on scales going non-linear (more so with biasing), 
it is not hard to get it more or less correct with ZA or ZA-based approximations. 

Coles, Melott and Shandarin (1993) (hereafter CMS) and Melott, Pellman and Shan- 
darin (1994a) (hereafter MPS) conducted a series of tests of analytical approximations by 
cross-correlating them with N-body simulations. They tested the Eulerian linear pertur- 
bation solution and ZA, i.e., the Lagrangian linear solution. They found that the latter is 
considerably more successful than the former. Considerable further improvement is made 
if ZA is applied not to the full power spectrum, but only to a truncated body of the spec- 
trum (TZ A) . The truncation removes unwanted non-linearities, which are evolved beyond 
the range of applicability of the model. The shape of high-frequency filters imposed on the 
power-spectrum was found to be important (MPS). Most successful in respect of improve- 
ment of the cross-correlation statistics was the use of Gaussian filters, but if one wishes 
to work backwards from evolved to initial state, a sharp truncation (step function) in k is 
preferable (Melott 1993). 

Most recent work has applied these tests to compare TZA with the "adhesion model" 
and the "frozen-flow-approximation" . The result was that those procedures provide rea- 
sonable statistical toy models, however, they are both dynamically less accurate than TZA 
as demonstrated by the comparison of their cross-correlation statistics in the fully devel- 
oped non-linear regime (Melott et al. 1994b, c). The comparison of higher moments - skew- 
ness and kurtosis - of the density contrast in the weakly non-linear regime (Bernardeau 
1994, Munshi & Starobinsky 1994, Munshi et al. 1994) has also demonstrated that the 
"frozen-flow-" and "frozen-potential-" approximations are dynamically less accurate than 
ZA. However, the results of the latter tests cannot be extrapolated to late non-linear stages, 
where also ZA has to be replaced by TZA. 

Because of the success of TZA it is worthwhile to study higher-order corrections to 
ZA as a subclass of the first-order Lagrangian perturbation solution (Buchert 1992) and 
also truncate the initial data in these models to see whether the performance of TZA can 
be further improved. We here continue our study on the Lagrangian perturbation schemes 
begun by Buchert, Melott and Weifi (1994a), henceforth BMW. In that paper we concen- 
trated on pancake models, i.e., models which do not involve small-scale power in the initial 
conditions. We did this as a first step in understanding the limitations of higher-order 
corrections to Zel'dovich's mapping. Also, pancake models can be understood as generic 
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archetype of hierarchical models which involve pancaking on all spatial scales (Kofman et 
al. 1992). In a second step, we here analyze (in the spirit of MPS) various power-spectra 
which are evolved deeply into the non-linear regime in order to probe the performance 
of the Lagrangian approximations in the case of hierarchical models. Although, it is not 
expected that a perturbation approach can be used to model highly non-linear stages, 
we here demonstrate that one can obtain good agreement with the N-body simulation by 
smoothing the initial data before evolving them. 

BMW found that the second-order approximation provides in pancake models a def- 
inite and useful improvement upon the first-order approximation (ZA). We explore now 
whether this improvement carries over to hierarchical models. 
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2. Numerical realization, Lagrangian perturbation schemes 
and cross-correlation statistics 



2.1. Numerical realization 

We specify initial data in terms of the power-spectrum V(k) (as a function of comoving 
wavenumber k = \k\) of Gaussian density perturbations of the form: 

V{k) =< \5 % \ 2 >oc \k\ n , (1) 

where 5^ is the discrete spatial Fourier transform of the density contrast 5, and < ... > 
denotes ensemble average over the entire distribution, assumed to be equivalent to the 
spatial average. We consider the range of indices — 3 < n < +1; for all spectra we use the 
same set of random phases, so that the overall structure is similar in all realizations. 

We emphasize that we have to give initial data early in order to fairly estimate the 
effect of the higher-order schemes and to guarantee an objective modelling of the collapse 
of first objects in the model (see the discussion in BMW). For the spectra with indices 
n = —3, —1 and n = +1 we give initial data at the r.m.s. density contrast of o~i(kN y ) = -05), 
where k^ y is the Nyquist frequency. For the other two indices we have taken a larger 
amplitude (cr^/c^) = 0.2), which is about that amplitude people normally use. This we 
did in order to see, whether there is any detectable effect. Indeed we found an effect which 
will be reported in Section 3 (see also Melott et al. 1994c). 

In order to let the concrete normalization of the spectra open, we define the evolution 
stages of the realizations in terms of the non-linearity scale k n i(t): 

rk n i(t) 

a 2 (t) / d 3 kV(k) = 1 , (2) 
Jo 

where k n i(t) is decreasing with time as successively larger scales enter the non-linear 
regime; a(t) is the scale factor of the homogeneous background (a(ti) := 1). The evolution 
of inhomogeneities is modelled in a flat Einstein-de Sitter background (a(t) = (j-) 2 ^ 3 ). 
We shall study two stages corresponding to k n i = 8k f and k n i = 4k /, where kf is the 
fundamental mode of the simulation box. Thus, our models studied here have been evolved 
for expansion factors of 240 to 5100, depending on spectral index and k n i. We shall present 
all statistical results for both stages. 

We evolve the fields with an enhanced PM (Particle-Mesh) method (Melott 1986) 
using 128 3 particles each on a comoving 128 3 mesh with periodic boundary conditions. 
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This method makes the code resolution-equivalent to a traditional PM code with 128 3 
particles on a 256 3 mesh (see Park 1990, 1991, and Weinberg et al. 1994). We then bin 
the data into a 64 3 mesh using the CIC scheme. Due to our controls (the same initial 
phases in all simulations) it is not necessary to do an ensemble of simulations with each 
power spectrum. One of each spectral type is sufficient to uncover systematic effects due 
to changing approximations. For more details see CMS. 
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2.2. Lagrangian perturbation schemes 

For the Lagrangian perturbation schemes up to the third order and their realization see 
Buchert (1994) and BMW. However, we here restrict the presentation of results to the 
first- and second-order schemes, the second-order scheme is expected to contain the major 
effects on large scales according to the results of BMW for the present purpose. But, we 
also made a consistency test by going to the third order. 

— * 

Denoting comoving Eulerian coordinates by q and Lagrangian coordinates by X, the 
field of trajectories q = F(X,t) up to the second order reads: 

F = X + qi (a)V S^(X) + q 2 (a)V S^(X) , (3) 

with the time-dependent coefficients expressed in terms of the expansion function a(t) = 
(f) 2:i : 




The perturbation potentials have to be constructed by solving iteratively the 2 Poisson 
equations: 

A <S (1) = I(S,i, k )ti , (3c) 
A SW = 2/1(5$ ) , (3d) 

where I and II denote the first and second principal scalar invariants of the tensor gradient 

I(S%)=tr(sV) = A SM , ( 3e ) 

^S) = ^[(^(^)) 2 -^((^) 2 )] • (3/) 

The general set of initial data is restricted according to the assumption of parallelism of the 
peculiar-velocity field and the peculiar-acceleration field at the initial time ti (see Buchert 
1994 and ref. therein for a discussion of this restriction). Therefore, the initial fluctuation 
field can be specified to the peculiar-velocity potential S alone. We can set = Sti 
(which is the unique solution of the first Poisson equation (3c) for periodic initial data, 
see Buchert 1992, Ehlers & Buchert 1994). Doing this, the flow-field (3) reduces to ZA if 
restricted to the first order. 
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We realize the solution by first solving Poisson's equation for S via FFT (Fast Fourier 
Transform) from the initial density contrast 6 generated as initial condition for the numer- 
ical simulation. In an Einstein-de Sitter model we have: 

AoS = ~6. (4) 

We then calculate the second principal invariant II directly from S and solve the second 
Poisson equation (3d) using FFT. The density in the analytical models is calculated by 
collecting trajectories of the Lagrangian perturbation solutions at the different orders into 
a 64 3 pixel grid with the same method (CIC binning) as in the N-body simulation. 

The CPU times on a CRAY YMP are for the first-order scheme 25 seconds, and for 
the second-order scheme 60 seconds; the corresponding CPU times on a CONVEX C220 
are 2 and 5 minutes. Thus, even the second-order scheme is competitive with one step in 
a corresponding PM-type N-body simulation. 
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2.3. Cross-correlation statistics 

For the presentation of the cross-correlation statistics we refer the reader to CMS, MPS 
and BMW for details. We use four statistics. Firstly, the usual cross-correlation coefficient 
S to compare the resulting density fields: 

S:=iMi, (5) 

(71(72 

where 5g,£ = 1,2 represent the density contrasts in the analytical and the numerical 
approximations, respectively, ag = a/ < Sg > — < Sg > 2 is the standard deviation in a 
Gaussian random field; averages < ... > are taken over the entire distribution. We believe 
this is the most important statistical test, because it measures whether the approximation 
is moving mass to the right place, with an emphasis on dense regions. We allow for small 
errors by presenting S for the two density arrays smoothed at a variety of smoothing 
lengths. 

Secondly, the power spectrum (eq. (1)) of the evolved N-body model and the analytic 
approximations were calculated. Thirdly, the mass density distribution for both of them 
are also shown. Both of these are widely displayed diagnostics in cosmology. 

MPS found that the primary reason for the superiority of one window function over 
another (equations (6) below) lay not in the final power spectral amplitudes but in the 
phase angle accuracy. We therefore follow MPS and will display < cos# where 9 = 
(pi — (p2 is the difference in the phase angle of the Fourier coefficients of mass density 
between the approximation and the simulation. The averaging is over all coefficients with 
the same magnitude of wave vector, as in the power-spectrum. 
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3. Filtering the initial data 



For the analytical realizations we filter the power-spectrum on various scales k c using 
a set of "trial windows". The use of Gaussian windows proved to yield the best cross- 
correlation with the N-body simulations in previous work on TZA (see MPS). We here 
again test three main types, sharp k— truncation (a), Gaussian windows (b), and tophat- 
smoothing (c). The filters (a) and (b) are described in Fourier space, whereas the filter (c) 
is applied in real space (but the computing for (c) is done in Fourier space). The initial 
data are modified to 5* = W5r, where the different filters W are: 

k K' 

W,Ak;k tr ) = {l : (6a) 

W G (k; ka) = e~ k 'l^ (66) 
... ,, _ , ( sm(R, h k) cos(R th k) \ 

Wth(hRth) = {lR^--(R^F) ■ (6c) 

We explored a wide range of values for the scale of each possible filter to find the optimum 
scale k op t for each type, as done for first-order (ZA) by MPS. 

In Figs.l we display the cross-correlation coefficient (eq. (5)) as obtained from a com- 
parison of the truncated analytical schemes with the raw data of the numerical simulations 
against the truncation scale for all spectra and for all filter types. This shows how the 
optimal truncation (defined as the scale with the largest cross-correlation coefficient) is 
varying with scale. In the cases with negative power index we notice that the variation in 
S with scale is much smaller towards higher k than in the cases with positive power index. 
The variation effect is very weak in cases with small-scale power index less than n = — 1. 

While type (c) gives much better agreement with the N-body results than simple 
k— truncation (a), we find as in MPS that the best filter type is always (b). Although 
the values for S for the best truncation are numerically similar for type (b) and (c), the 
truncation scale for the best agreement is always larger for type (c), thus erasing more 
non-linear information. We therefore concentrate our further analysis on type (b) and 
recommend generic use of Gaussian filtering. 

In Table 1 we display the values of S for the optimal truncation scales k opt for the 
different power-spectra and Lagrangian schemes. In Fig. 2 we illustrate this table by plot- 
ting S versus the power index n. We see an almost constant (but small) improvement of 
the second-order scheme over first-order for positive indices. As soon as the index drops 
negative we find increasing improvement with decreasing power. Fig. 2 also shows that 
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we obtained slightly different performance for the models with spectral index n = and 
n = —2. These models have been started at a later stage (compare 2.1). We think that this 
could be an indication for inaccuracies of the numerical evolution code for low-amplitude 
initial conditions, which we want to test in future work (see also Melott et al. 1994c). 

We have used discrete steps in k to find the optimal truncation lengths. If we express 
the optimal scales in terms of the non-linearity scale, e.g., k n \ = 8/c/, we find that the 
arithmetic mean (used here as a spectrum independent estimate) lies for the first-order 
scheme at 1.45fc n /, for second-order at l.2k n i, i.e., at a slightly larger scale. This difference 
almost vanishes at the later stage (k n i = 4/c/), where the arithmetic means are 1.5/c n ; for 
first-order and lAk n i for second-order. Still, within the uncertainty bounds involved in 
the determination of the optimal scale we can recommend generic use of a k— scale close 
to, but larger than the non-linearity scale as suggested by MPS. In real space this means 
that the comoving truncation scale At < X n i- It is interesting that, although the second- 
order scheme requires slightly larger truncation scales, the dependence of that scale on the 
spectrum is weaker than in the first-order scheme. We emphasize that for negative-sloped 
spectra the truncation scale At is especially small; in view of Fig.l a truncation of models 
with negative spectral slope is not as important as in the other models. For spectra n < —2 
the differences in S on the small-scale end are negligible. These models perform well even 
without truncation as will become evident from the Figs. 4 below. 

We proceed by statistically analyzing the optimally truncated models only, with k opt ac- 
cording to Table 1. The results are presented and discussed in the next section. 
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4. Discussion of the results 



We first concentrate on the first stage corresponding to k n i = 8k f. 

Fig. 3 presents the cross-correlation coefficient with the N-body simulation as a func- 
tion of scale for the optimally truncated and untruncated first- and second-order schemes 
together with two versions of the (still widely used) Eulerian linear approximation: (1) 
linear density contrast, and (2) (following CMS): "chopped" linear density contrast, i.e.: 
1 + ^chopped = ot(l + 5), if 5 > — 1 and otherwise, where a is the normalization constant 
keeping the total mass the same. We have chosen the spectra with indices n = — 1 and 
n = +1 for the illustration of the differences among the different approximations. This 
comparison again demonstrates that the Lagrangian schemes (even in their untruncated 
form) perform substantially better than the Eulerian solution with or without "chopping" . 
In comparison with the following figures this demonstration should teach that the Eulerian 
approach is almost useless to describe even slightly non-linear ly evolved density fields. 

In Figs. 4 we display slices of the density contrast field as predicted by the N-body 
simulations and the two optimally truncated Lagrangian perturbation schemes. We also 
display the evolved fields for the untruncated perturbation schemes to illustrate the effect of 
truncation, and for the "chopped" Eulerian first-order model. The structures for different 
power-spectra bear family resemblance due to our phase controls. 

Generally, we appreciate a reasonable correspondence of the numerical density fields 
with those of the truncated Lagrangian schemes. Although, in the cases with small amount 
of power on small scales (especially n = —2, —3) the pictures for the untruncated models 
might suggest better coincidence with the numerical density fields, the cross-correlation 
statistics with the raw numerical data (shown in Figs.l) prove that the truncated models 
agree slightly better. However, in these models the truncation plays a peculiar role, since 
the cross-correlation coefficient only varies weakly as a function of scale by going to smaller 
scales. The small numerical differences in the coefficients can be seen in Figs.l. Already 
small differences in the N-body computing can yield different "optimal" truncation scales. 
Within this uncertainty, it is a matter of standards we want to apply, whether we truncate 
the models n = —2, 3 or not. But, clearly for models with spectra n > —2, truncation of 
the initial data is indispensible. 

The contours of the optimally truncated first- and second-order models are quite 
similar; visually manifest differences only occur for the cases n = — 2 and n = —3. Note 
that this is also due to the presentation in terms of grey-scale plots which are truncated 
at small density levels (density contrast 5), where both schemes predict similar contours. 
The higher cross-correlation coefficient of the second-order scheme reported below is the 
result of sharper and higher density peaks at second order, visible at higher density levels. 
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In Figs. 5 we show the r.m.s. density contrast of the N-body simulations together with 
that of the analytical schemes as a function of scale. In the following figures we shall use 
the r.m.s. density contrast of the N-body simulations as a scale parameter. 

Figs. 6 show the most important result of this paper: the cross-correlation coefficient 
for different smoothing scales (smoothed with a Gaussian window) for different power- 
spectra and the two truncated Lagrangian schemes. We infer that the second-order ap- 
proximation shows higher cross-correlation with the N-body density field for all spectra 
and for all scales. However, at this stage the improvement is small in cases with non- 
negative power index, but it grows substantially with decreasing small-scale power. This 
result can already be seen in Fig. 2. In general, we find that for negative power indices 
the agreement with the N-body field is very good, especially on scales which are rele- 
vant for modelling large-scale structure (i.e., above the scale of galaxy clusters), bearing 
in mind that both the numerical simulations and the analytical schemes are assumed to 
be approximations to the unknown exact solution of the problem. On scales far below 
the non-linearity scale, where the r.m.s. density contrast of the N-body simulation as- 
sumes the value 2, the cross-correlation coefficient generally lies close to S = 0.8 in the 
first-order scheme, while the second-order scheme is closer to S = 0.9 for negative-sloped 
power-spectra. 

While the cross-correlation test yields information on the dynamical capabilities of 
the models, the following statistics show whether they are useful in reproducing statistical 
properties of the density field. 

In Figs. 7 we display the power-spectra as calculated from the final distributions. 
We confirm the results by MPS who found a weak representation of the power-spectrum 
on scales below the non-linearity scale in TZA for positive power-spectra slopes. Here, 
the second-order result does not improve upon first-order. Again, for negative slopes, the 
agreement is reasonably good for both approximation schemes also below the non-linearity 
scale. In summary, for much small-scale power the power-spectrum is weakly represented 
below the non-linearity scale and especially below the co moving truncation scale At- This 
is evident, since a truncation of the spectrum removes power which cannot be compensated 
by non-linear evolution. 

The power-spectrum provides low-order statistical information. Figs. 8 and 9 present 
statistics which probe higher-order information of the density distribution. We find that 
phase-information is represented much better in the second-order scheme in the cases 
n = — 2 — 3, while in all other cases the performance of second-order is moderately better 
than the first-order model. Also here we see that the positive second-order effect is 
counteracted by a large amount of small-scale power at this stage which erases phase- 
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information present in the weakly non-linear regime. 

The density distribution functions (Figs. 9) show excellent agreement with the N-body 
result for all cases with negative power index; for n = 0, +1 the number of high-density cells 
is underestimated and the number of low-density cells is overestimated by both schemes. 
The improvement due to the second-order correction is still measurable. 

For the sake of showing consistency of our results, we made two further studies: (1) 
we evolved the fields to a later non-linear stage, where the non-linearity wavenumber has 
dropped to half the value of the stage analyzed above; (2) we also calculated the statistics 
for the third-order model (see BMW). 

In Figs. 10-14 the statistics corresponding to Figs. 5-9 are displayed for the later stage 
with k n i = 4k f. 

One of the striking results obtained by MPS concerning the cross-correlation analysis 
is that TZA has proved to be very robust, i.e., the good performance of this model remains 
stable by going to later stages, or to smaller scales, respectively. In the context of the 
Lagrangian perturbation solutions, BMW attributed this property to the fact that TZA is 
the principal first part of a perturbation sequence, while higher-order corrections display a 
"blow-up" phenomenon at and after the stage when the perturbation theory breaks down. 
Interestingly, we here find that the improvement found for the second-order Lagrangian 
scheme is also robust, i.e., at later stages and on small scales the truncated second-order 
model does not become worse than TZA as observed for pancake models by BMW. This 
feature can be traced back to the optimal filtering employed here which compensates the 
"blow-up" of the higher-order coefficients by the choice of a different optimal scale for each 
order. Moreover, the improvement of second- upon first-order is even larger at the later 
stage, especially with respect to the cross-correlation coefficients, the phase-errors, and the 
density distribution functions. In particular, also the spectra with positive power-index 
show improvement, especially on small scales. 

Since from previous comparisons with other approximations we know that it is hard 
to improve on the dynamical performance of TZA, we consider this as a definite and 
useful success of the second-order model. It is remarkable that, although the second-order 
scheme needs truncation at a slightly smaller frequency than first-order, thus erasing more 
information initially (compare Fig.l), the non-linear evolution overcompensates that. 

We could not confirm the same property of 'robustness' for the third-order scheme. 
In BMW we already noted that differences between second- and third-order are smaller 
than the difference found for second- upon first-order, a result which is expected for 
a perturbation approach which converges to the exact solution (if it does). However, 
the break-down phenomenon at and after the perturbation schemes are evolved beyond 
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the range of validity after shell-crossing, i.e. the behavior observed in pancake models 
by BMW, could be compensated by optimized filtering of initial data for second-order 
perturbations. This seems not to work for the third-order scheme. Here, we have to stress 
that the numerical realization of the third-order scheme by iteratively solving Poisson 
equations should be very accurate in order to draw definite conclusions, especially for 
models with much small-scale power. We think that much more numerical effort has to 
be added in order to treat the third-order effect reliably. Remember that in the third- 
order model 7 Poisson equations have to be solved, 4 of them based on solutions of the 
second-order Poisson equation, while in the second-order model only 2 Poisson equations 
have to be solved (directly from the initial data). Buchert et al. (1994b) therefore pursue 
a different way to overcome numerical problems by solving Poisson equations analytically 
up to the third order. With this method one can obtain analytical expressions for all 7 
perturbation potentials, and realize the Lagrangian scheme for random fields with up to 
about 100 Fourier modes at reasonable CPU speed. First results of this comparison are 
discussed by Buchert (1993b). 

We emphasize that a breakdown of higher than second-order solutions at late non- 
linear stages does not contradict the result that third-order clearly improves upon second- 
order in the weakly non-linear regime before shell-crossing, e.g. in comparison with the 
exact spherically symmetric solution, or the skewness and kurtosis moments of the density 
field (Munshi et al. 1994); here we tested the extreme situation of going to stages which 
involve a large number of shell-crossings. 

We finally remark that that we were using the exact expressions for second- and 
third-order solutions for initial conditions, where the peculiar- velocity potential S and 
the peculiar-gravitational potential are initially proportional (S = —<j>ti). The expres- 
sions used by other authors on higher-order Lagrangian schemes refer to different initial 
conditions and do not involve terms other than the leading term. The solutions we use 
involve growing parts (even with different signs than the leading terms) due to second- and 
third-order homogeneous solutions of the basic differential equations. The leading terms 
are the particular solutions of these equations only. Of course. ctt> lctte non-linear stages an- 
alyzed here, only the leading terms matter, which we confirmed by dropping all terms other 
than the leading terms. We found that for expansion factors a < 25 differences with the 
calculations presented here become significant (compare, e.g., the second-order coefficient 
of the particular solution of the displacement field — 3/14a 2 with the growing homogeneous 
solution +3/5a; at a ~ 25 the difference still is a 10% effect in the time-coefficients). 
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5. Conclusions 



We have analyzed a family of hierarchical models with different amount of power on small 
scales, and evolved deeply into the non-linear regime of structure formation. We have 
compared the results from N-body simulations with those from Lagrangian perturbation 
schemes, where for the latter we have filtered the initial data before evolving them. The 
filter scale was determined by optimization of the cross-correlation coefficient with the 
corresponding N-body simulation. 

We can report two striking advantages of going to second order in perturbation the- 
ory. Firstly, especially the statistics which probe the gravitational dynamics of the models 
show improvement due to second-order corrections. This success is found for a consider- 
ably higher non-linearity than is expected from a perturbation approach. Secondly, the 
improvement (although small for much small-scale power) is robustbj going to later stages 
and to smaller scales. 

We conclude that the second-order scheme with truncated power-spectra will be 
especially useful for modelling large-scale structure. In particular, the second-order scheme 
is computationally simple to realize with a CPU speed comparable to that of TZ A. Since the 
second-order corrections to TZA provide noticeable improvement of dynamical accuracy 
for initial data with negative sloped power-spectra, we expect that the truncated second- 
order scheme will be useful for the modelling of standard cosmogonies (like Hot-, Cold-, 
and Mixed-Dark-Matter). This modelling will be especially effective for large sample 
calculations, since in numerical realizations of 'fair' samples in excess of 300 h _1 Mpc, 
performed with the same resolution as the simulations reported here, the truncation scale 
is close to the Nyquist frequency of the N-body computing. Thus, shortcomings of the 
analytical schemes become negligible which puts them in an ideal position for the purpose 
of simulating the environment of galaxy formation down to scales, where other physical 
effects start to affect models which are based on the description of self-gravity alone. Our 
method can be effective down to galaxy group mass scales (10 13 M Q ), or better if we include 
biasing or go to epochs earlier than the present. 
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Figure Captions 



Figure 1: The variation of the cross-correlation coefficient S (eq. (5)) against the truncation 
scale k c imposed on the initial data for the evolved analytical schemes (k n i = 8k f) for different 
filter types (sharp k— truncation (kc; eq. (6a)), Gaussian truncation (gs; eq. (6b)), and (real space) 
top-hat truncation (th; eq. (6c))), for spectral indices n = + 1, 0, — 1, — 2, — 3 (Fig.la,b,c,d,e). 
The first-order model is shown as a dotted line, the second-order model as a dashed line. 

Table 1: The cross-correlation coefficient S for the optimal truncation scale k op t in units of 
k n i = 8k f (stage f), and k n i = 4k f (stage g) for different spectral indices and the two Lagrangian 
perturbation schemes. 

Figure 2: Illustration to Table 1. Shown is the cross-correlation coefficient S for the optimal 
truncation scale k op t according to Table 1 for stage f (left) and stage g (right) as a function 
of spectral index. The first-order scheme is shown with squares, the second-order scheme with 
diamonds. 

Figure 3: The cross-correlation coefficient S for the stage k n \ = 8k f and the cases n = — 1 
(upper) and n = +1 (lower) as a function of scale. Displayed is the cross-correlation of the 
N-body simulation with the Eulerian linear extrapolation of the initial data (dashed-dotted) , the 
"chopped" Eulerian linear approximation (thick dashed-dotted), the untruncated first- (dotted) 
and second-order (dashed) approximations, and the optimally truncated first- (thick dotted) and 
second-order (thick dashed) approximations. 

Figure 4: A thin slice (thickness L/64) of the density field is shown for the numerical (upper- 
left), the "chopped" Eulerian linear (upper-right), the optimally truncated first-order (middle- 
left) and second-order (lower-left), the untruncated first-order (middle-right), and the untrun- 
cated second-order (lower-right) approximations for the evolution stage corresponding to k n i = 
8kf, and for different power-spectra with index n = +1, 0, — 1, —2, —3 (Figs. 4a,b,c,d,e). The 
grey-scale is linear, and the maximum density contrasts are adapted to give a satisfactory visual 
impression, except for the optimally truncated models for which the upper truncation is the same 
(density contrast 5) as for the numerical simulation to allow for an objective comparison. 

Figure 5: The standard deviations of the density contrast as a function of smoothing scale 
(given in grid units Rq of a 64 3 grid) in the optimally truncated first-order (dotted), the optimally 
truncated second-order (dashed), and the numerical (full line) approximations, for the different 
power-spectra n = —3, —2, —1, 0, +1 (Figs. 5a,b,c,d,e,). 
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Figure 6: The cross-correlation coefficient S as a function of the standard deviation a p of the 
numerical simulation for the different power-spectra n = — 3, — 2, — 1, 0, +1 (Figs. 6a,b,c,d,e). 
The cross-correlation of the N-body with the optimally truncated first-order model is shown as 
a dotted line; with the optimally truncated second-order model a dashed line. 

Figure 7: The power-spectra of the N-body simulations (solid line) compared with the opti- 
mally truncated first- (dotted) and second-order (dashed) approximation schemes for the different 
power-spectra n = —3, —2, — 1, 0, +1 (Figs. 7a,b,c,d,e). 

Figure 8: The relative phase-errors. Notation and labelling like in Figure 6. 

Figure 9: The density distribution functions. Notation and labelling like in Figure 7. 

Figure 10: Same as Fig. 5 for a later stage when the non-linearity scale has dropped to half the 
value of the former stage (k n i = Ak f). 

Figure 11: Same as Fig.6 for the later stage. 

Figure 12: Same as Fig. 7 for the later stage. 

Figure 13: Same as Fig.8 for the later stage. 

Figure 14: Same as Fig.9 for the later stage. 
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